Abstract. We present detailed thermal and gas-phase 
chemical models for the envelope of the massive star- 
forming region AFGL 2591. By considering both time- and 
space-dependent chemistry, these models are used to study 
both the physical structure proposed by van der Tak et al. 
( |1999| ; |2000| ), as well as the chemical evolution of this re- 
gion. The model predictions are compared with observed 
abundances and column densities for 29 species. The ob- 
servational data cover a wide range of of physical condi- 
£N| tions within the source, but significantly probe the inner 
. regions where interesting high-temperature chemistry may 
' be occurring. Taking appropriate care when comparing 
" , models with both emission and absorption measurements, 
we find that the majority of the chemical structure can 
be well-explained. In particular, we find that the nitrogen 



and hydrocarbon chemistry can be significantly affected 
[■ — , by temperature, with the possibility of high-temperature 
t-H ' pathways to HCN. While we cannot determine the sul- 
phur reservoir, the observations can be explained by mod- 
' els with the majority of the sulphur in CS in the cold gas, 
, SO2 in the warm gas, and atomic sulphur in the warmest 
' gas. Because the model overpredicts CO2 by a factor of 40, 
f\J , various high-temperature destruction mechanisms are ex- 
' plored, including impulsive heating events. The observed 
, abundances of ions such as HCO+ and N 2 H+ and the 
' cold gas-phase production of HCN constrain the cosmic- 
--^ , ray ionization rate to ~ 5.6 x 10~ 17 s _1 , to within a factor 
1-^ ' of three. Finally, we find that the model and observations 
can simultaneously agree at a reasonable level and often to 
O ' within a factor of three for 7 x 10 3 < t(yrs) < 5 x 10 4 , with 
+3 , a strong preference for t~3x 10 4 yrs since the collapse 
^ ' and formation of the central luminosity source. 



Key words: Stars: formation - Stars: individual: AFGL 
2591 - ISM: molecules 



A&A manuscript no. 

(will be inserted by hand later) 

Your thesaurus codes are: 

missing; you have not inserted them 



Chemistry as a probe of the structures and evolution of 
massive star-forming regions 

S. D. Doty 1 , E. F. van Dishoeck 2 , F. F. S. van der Tak 2 - 3 , and A. M. S. Boonman 2 

1 Department of Physics and Astronomy, Denison University, Granville, OH 43023 

2 Sterrewacht Leiden, P.O. Box 9513, 2300 RA Leiden, The Netherlands 

3 Max-Planck-Institut far Radioastronomie, Auf dem Hiigel 69, 53121 Bonn, Germany 

Received; accepted 



ASTRONOMY 

AND 
ASTROPHYSICS 



1. Introduction 

The distribution and composition of dust and gas around 
isolated low-mass young stellar objects (YSOs) has been 
well-studied both observationally and theoretically. Unfor- 
tunately, much less is known about the distribution and 
composition of material around high-mass YSOs (see e.g., 
Churchwell 1993, 199E). The higher densities and masses, 
and shorter lifetimes associated with massive star forma- 
tion suggest that differences between regions of high- and 
low-mass star formation can be expected. 

Recent observational advances (e.g., submillimeter 
beams of ~ 15" sampling smaller regions of higher critical 
densities, interferometry at 1 and 3 mm, and ground- and 
space-based infrared observations of gas and ices) have 
led to a new and better understanding of the environ- 
ment around massive YSOs (see e.g., Garay & Lizano 



199S, van Dishoeck fc Hogerheijde 1999, Hatchell et al. 



200C, Beuther et al. 2002). In this vein, van der Tak et 



al. ( 199£ , 2000 ) have conducted detailed multi- wavelength 
studies of high-mass YSOs, and begun to form a picture 
for the physical structure of some of these regions. The 
proposed material distributions in the envelopes fit a wide 
variety of continuum and spectral line data. However, they 
are incomplete without a detailed thermal and chemical 
structure. The proposed material distribution can be used 
to test the chemical structure and evolution of the enve- 
lope, and the combined results can eventually be used to 
compute line strengths and profiles for direct comparison 
with observations. 

Significant work has been involved in developing an un- 
derstanding of the chemistry of star- forming regions. This 
ranges from studies of cold, dark clouds (e.g., Herbst & 
Klemperer 1973, Prasad & Huntress 1980, Leung, Herbst, 



& Huebner 1984, Gwenlan et al. 200C) to "hot cores" (see 
e.g., reviews by Millar \1993\ Walmsley & Schilke |1992; 



Kurtz et al. 2000). In nearly all cases, however, the chem- 
istry is considered for a homogeneous cloud, or a point 
within a cloud (see, though, Xie et al. 1995 and Bergin et 
al. 1995| for counter examples). Unfortunately, the physical 



conditions (i.e., temperature and density) vary strongly 
with position within the envelope, meaning that poten- 
tially extreme chemical variations may occur between the 
source center and the observer. It is this strong variation 
of chemical composition with position and time that may 
provide one of the best benchmarks of our understand- 
ing of both the structures and evolution of massive star- 
forming regions. 

In this paper, we utilize position-dependent thermal 
balance and time- and position-dependent chemical mod- 
eling to probe the validity of the physical structures pro- 
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posed by van der Tak et al. (1999 & 2000), and more im- 
portantly, to study the chemical evolution of AFGL 2591. 
In particular, taking their structure as a starting point, 
we construct detailed models for the gas-phase chemistry 
of this source, and compare the results with observations. 
AFGL 2591 is a massive (~ 42M within r = 3 x 10 4 AU), 
luminous (~2x 10 L@) infrared source with many of the 
properties thought to characterize YSOs. While most mas- 
sive stars form in clusters, AFGL 2591 has the advantage 
that it is forming in relative isolation - allowing us to 
study its physical, thermal, and chemical structures with- 
out influence from other nearby massive sources. It has 
the further advantage of being well-observed both in the 
continuum and in a variety of molecular lines. 

This paper is organized as follows. The existing ob- 
servations providing the model constraints are briefly dis- 
cussed in Section 2. In Section 3, the model is described. 
The model is then applied to AFGL 2591 and compared 
with the observational results in Section 4. In Section 5, we 
compare our time-dependent model predictions with the 
observations in order to constrain the chemical age of the 
envelope. Finally, we summarize our results and conclude 
in Section 6. 



2. Existing observations and usage 

AFGL 2591 has been well-observed both in the continuum 
and in various molecular lines. While no new observations 
are presented in this paper, it is important to briefly note 
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and discuss the relevant observations as they provide the 
constraints placed on the model. 

2.1. Continuum 

AFGL 2591 has been observed in the range 2 — 60000 fxm 
by Lada et al fll984| ), Aitken ( |198SD , Sand ell (19 98, private 
communication), and van der Tak et al. ( 1999[) . These re- 
sults were analyzed by van der Tak et al. ( ]2000| - see Sect. 
3 below) to constrain the density distribution and grain 
properties - necessary for not only the thermal structure, 
but also to properly evaluate the gas thermal balance and 
hence obtain the gas temperature as a function of position. 

2.2. Molecular lines 

A wide variety of observations, both in the infrared and 
submillimeter, have been conducted of molecular gas in 
AFGL 2591, some of which are as of yet unpublished. The 
results are summarized in Table |l|, where the species, ob- 
served abundance [x(X) = n(X)/n(H 2 )] or column density 
[Af(X)], inferred excitation temperature, method of anal- 
ysis, weight used in selecting the most important of the 
relevant observations, type of observation, and reference 
are listed. 

2. 3. Notes on Table [I] and usage of the data 

The observation type is listed in Table [l] as this is 
significant for comparing the results with observations. 
For infrared absorption lines, the molecules observed are 
along the (narrow) line of sight to the background con- 
tinuum source. Consequently, these results should be 
compared to model "radial column densities", namely 
-^radial = Jn(r)dr. On the other hand, submillimeter 
emission lines arise from throughout the envelope. In 
these cases, averages over the beam are used in compar- 
ing predicted and observed column densities. Here, the 
"beam-averaged column density" is defined as JVb cam = 
/ / n(z,p) dz G(p)2irp dp/ J G(p)2irp dp where p is the 
impact parameter, and G(p) is the beam response func- 
tion. We also divide the data in this fashion, as we expect 
many of the uncertainties in the analysis to be similar for 
one type of observation. 

In columns 2 and 3 of Table [l] we list the inferred frac- 
tional abundance or column density of the given molecule 
toward AFGL 2591. This is done to provide the most com- 
prehensive set of information with which to compare our 
models. 

While determination of column density is relatively 
straightforward for infrared absorption lines in the limit of 
no re-emission, the situation is more difficult for emission 
lines as the emission may arise from a range of radii, and 
thus a range of densities, temperatures, chemical abun- 
dances, and optical depths. To combat this, some ef- 
fort has been made recently to determine the fractional 



abundance within the envelope through detailed, non lo- 
cal thermodynamic equilibrium (NLTE) radiative trans- 
fer (RT) modeling (van der Tak et al. |1999[ ). When this 
is done, we view the inferred abundances as superior to 
pure column densities as they account for many of the 
potential errors in determining the column density. As ex- 



amples, van der Tak, van Dishoeck, & Caselli (200C) and 
Boonman et al. ( |2001 ) have used such modeling to suggest 
"jump" models for the chemical enhancement of species 
within certain regions of YSO envelopes. As a result, in 
column 5 we note the method used in determining the 
observational result. We also use these criteria to assign 
a weight (higher is better) in column 6 to denote which 
data/fits we view as superior. In cases where the fit due to 
radiative transfer modeling is only moderate, we give this 
result the same weight as the results from other methods. 

In both cases, where an excitation temperature can be 
assigned to the data, we note the temperature for that 
component as T cx in column 4 of Table [|. While T ox is 
not necessarily equal to the kinetic temperature, it does 
give some indication as to the region from which the ob- 
servation arises. The values of 38, 200 and 1010 K refer to 
the excitation temperatures of CO found by Mitchell et 
al. (1989) in infrared absorption line studies. The 200 K 
component is thought to be associated with shocked out- 
flowing material, whereas the other two temperatures refer 
to the quiescent envelope. 

Finally, we note the relative importance of different 
measurements for probing various regions in the envelope. 
Absorption is confined to the narrow line of sight toward 
the central source. For centrally-condensed envelopes, the 
column density is dominated by the interior. This makes 
absorption measurements useful for probing the warm in- 
terior. On the other hand, emission measurements can and 
often do arise from throughout the envelope. When the 
density falls off slower than r -2 as is the case for AFGL 
2591 (van der Tak et al. 199S ) the outer portion of the en- 
velope dominates the mass, and so emission measurements 
are often more useful for probing the cool exterior. These 
expectations are relatively consistent with the results of 
Table |, where many of the absorption measurements in- 
clude significant high excitation temperature components, 
while the inferred excitation temperatures for the emission 
data are generally much lower. 



3. Model 

In this section, a brief synopsis of the physical, thermal, 
and chemical models are provided. For more d etailed in- 
formatio n, see van der Tak et al. (|l999j |2000|) , Doty & 
Neufeld ( 1997 ), and references therein. For reference, the 
model parameters are reproduced in Table ^[ 
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Table 1. Inferred Column Densities and Abundances Toward AFGL 2591 
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3.1. Physical model 

Our model for AFGL 2591 concentrates on the extended 
envelope of source. While an inner disk may be present, 
OVRO interferometric observations by van der Tak et 
al. ( 1999[ ) suggest an unresolved central source of radius 
30 < r(AU) < 1000. These and other continuum observa- 
tions were analyzed by van der Tak et al. (2000) using a 
modified version of the self-consistent continuum radia- 



tive transfer model of Egan, Leung, & Spagna ( |1988| ). 
Based upon the fit to the continuum flux and surface 
brightness, as well as CS line data, they constrained the 
density to a power law of the form n(r) = no(ro/r) a . 
In particular, they found a best fit with a = 1.0, and 
n = n(R 2 ,r = r = 2.7 x 10 4 AU) = 5.8 x 10 4 cm" 3 . We 



adopt these values for the remainder of the paper. Our 
inner radial position was chosen to be r- m = 2 x 10 2 AU, 
corresponding to T ~ 440 K (see below). This inner radius 
was chosen not only for consistency with the observations 
of van der Tak et al. ( |2000| ) , but also as extrapolation of a 
density power law further into the interior of the envelope 
leads to column densities inconsistent with observations 
(see Sect. 4.2 below). 



1999 



Following the analyses of van der Tak et al. 
pOOOD and Doty & Neufeld Ql997| ), we assume that the 



physical and thermal model does not change significantly 
with time so that an equilibrium may be achieved, but we 
do allow for a time-dependent chemical evolution. While 
the collapse and rise in luminosity will occur on short 
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time scales (< 1000 yr), it is the combination of density 
structure and luminosity of the central source that sets 
the temperature structure. Therefore, as long as the enve- 
lope mass and luminosity do not significantly change, we 
can consider the source as approximately constant over 
the ~ 10 5 yrs in which the envelope will be dissipated 



(Hollenbach et al. 1994, Richling & Yorke 1997). To see 



this, consider the fact that the free fall and sound-crossing 
times at the outer edge are both ~ 2x 10 5 yrs. While these 
timescales are smaller closer to the center, accretion events 
should probably only be important in the very interior. 

Finally, we note that an outflow has been observed to- 
ward this source (see, e.g., Bally & Lada 1983| , Mitchell 
et al. 1989). However, spectroscopy shows that nearly all 
submillimeter lines with the exception of CO can be as- 
signed to the envelope as their linewidths arc only ~ few 



km s _1 (see, e.g., van der Tak et al. 1999, and recent and 
upcoming infrared data from TEXES by Knez et al. 2002 
and Boonman et al. - in preparation). Only CO has a sig- 
nificant fraction of the observed material in the outflow. 
This assignment of material to the envelope rather than 
the outflow is also justified a posteriori as our models are 
able to reproduce a good deal of the observed chemistry 



without the requirement of shock chemistry. Because the 
submillimeter lines probe high excitation gas, the lower 
density surrounding cloud is automatically filtered out. 



3.2. Thermal model 

The equilibrium gas temperature within the cloud is de- 
termined by the balance between heating and cooling. The 
gas heating is dominated by gas-grain collisions, and the 
dust temperature is determined from the self-consistent 
solution to the continuum radiative transfer problem as 



above. The Neufeld, Lepp, & Melnick (1995) cooling func- 
tions were adopted, with modifications as noted in Doty 
& Neufeld (1997). Furthermore, as the Neufeld, Lepp, & 
Melnick (1995) cooling functions were constructed assum- 
ing a singular isothermal sphere (with a commensurate 
n oc r~ 2 density power law), they were modified to be 
applicable to the r _1 power law adopted here. This en- 
tailed two corrections. First, the column densities had to 
be computed correctly at each position, rather than sim- 
ply relying upon the local density. Second, the cooling 
functions for the tabular results were modified by a fac- 
tor [N(a = 2)/N(a)] f , where N(a) = n(r) dr is the 
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Table 2. Model Parameters 
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Fig. 1. Physical and thermal structure of AFGL 2591. The 
density and dust results from the model of van der Tak 
et al. ( 2000 ). The gas temperatures are calculated from 
t he de tailed thermal balance, similar to Doty & Neufeld 
( |l997|) . Note that T gas ~ T dust . 
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cm -2 , and / = — 1.0 for N > 10 21 cm" 2 . This factor is 
chosen to match the functional dependence of the cool- 
ing rate on the column density as described in Neufeld, 
Lepp, & Melnick ( 1995 ) for H2O - the dominant coolant 
utilized in tabular form - and is consistent with the fact 
that the cooling rate should be inversely proportional to 
the column density for opaque sources. The resulting gas 
temperature distribution is shown in Fig. S, and is physi- 
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cally similar to that of Doty & Neufeld (199 
T^as ~ Tdust, as was assumed by van der Tak et al 
200C ). For comparison, models run assuming T gas 



namely that 



( |1999 , 

= Tdust 



Parameter 


Value 


Ref. 


Outer radius (AU) 


3.0(4) 


a 


Inner radius (AU) 


2.0(2) 


a 


Density [n(r) = no(ro/r) a ] 






. . . Exponent [a] 


1.0 


a,b 


. . . Ref. position [ro] (AU) 


2.7(4) 


b 


. . . Ref. H2 density [no] (cm -3 ) 


5.8(4) 


b 


CR ionization rate [f] (s -1 ) 


5.6(-17) 


c 


Initial Abundance (T > 100K) 






CO 


3.7(-4) 


a 


CO2 


3.0(-5) 


d 


H 2 


1.5(-4) 


d 


X120 


1 f\f 

1 .u\ UJ 


oee texL 


N 2 


7.0(-5) 


e 


CH4 


1.0(-7) 


e 


C2H4 


8.0(-8) 


e 


C2H6 


1.0(-8) 


e 


OI 


0.0(0) 


e 


H2CO 


1.2(-7) 


c 


CH3OH 


1.0(-6) 


c 


S 


0.0(0) 


e 


Fe 


2.0(-8) 


e 


Initial Abundance (T < 100K) 






CO 


3.7(-4) 


a 


CO2 


0.0(0) 


f 


H 2 


0.0(0) 


f 


H 2 S 


0.0(0) 


f 


N 2 


7.0(-5) 


e 


CH4 


1.0(-7) 


e 


C2H4 


8.0(-8) 


e 


C2H6 


1.0(-8) 


e 


OI 


8.0(-5) 


g 


H 2 CO 


0.0(0) 


f 


CH3OH 


0.0(0) 


f 


s 


6.0(-9) 


see text 


Fe 


2.0(-8) 


e 



10" 



All abundances are gas-phase, and 



a(b) means a x 
relative to H2 
a van der Tak et al. 

Tak & van Dishoeck 2000, d Boonman et al. 



lOOq b van der Tak et a l. POOP 



200C 



c van der 
Charnley 

1997, f assumed frozen-out or absent in cold gas-phase, g taken 
to be ~ consistent with Meyer, Jura, & Cardelli 1998 



show no significant differences. 



3.3. Chemical model 

The chemical model is based upon the UMIST gas-phase 
chemical reaction network (Millar, Farquhar, & Willacy 
1997). Using this network, we construct pseudo time- 
dependent models of the evolution of the chemical abun- 
dances. We do this over a range of 30 radial grid points, 
providing a time- and space-dependent chemical evolution. 
The local parameters (hydrogen density, temperature, and 
optical depth) at each radial point are taken from the 
physical and thermal structure calculations above. For 
our initial abundances, we follow Charnley fll997| ; private 



communication). These parameters allow us to reproduce 
many of the results of the hot core models of Charnley 
(1997; private communication), with most discrepancies 
directly attributable to differences in adopted reaction 
rates. 

We also include the approximate effects of freeze-out 
onto dust grains by initially depleting certain species be- 
low 100 K (see Sect. 4.7 for discussion of H 2 CO and 
CH3OH). We attempt to minimize this effect by predom- 
inantly depleting those species that have high observed 
solid-phase abundances. Our initial fractional abundances 
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not in CO is initially placed into water. This is consistent 
with models we and others (e.g., Doty & Neufeld 1997, 



15.5 16 16.5 17 

Log[ Radial Position / cm ] 



17.5 



Fig. 2. The fractional abundances of CO and H2O 
throughout the envelope as a function of time. The 
dashed-lines correspond to the (constant) CO abundance, 
and the solid lines to the H2O abundance. The curves are 
labeled by the time in years, where a(b) = a x 10 b . 

relative to H2 , as well as other model parameters are listed 
in Table |. 

The cosmic-ray ionization rate is taken from van der 
Tak & van Dishoeck (|2000j ) for AFGL 2591, and will be 
discussed in Sect. 4.5. The effects of cosmic-ray induced 
photochemistry were ignored. The initial sulphur abun- 
dance was chosen to make the models agree with obser- 
vations (see Sect. 4.6). The assumed sulphur abundances 
are in general agreement with observations for both the 
warm (e.g., toward Orion by Minh et al. 
cold (e.g., Irvine, Ohishi, & Kaifu 1991) components. 



199C), and the 



The effects of photodissociation from the ISRF at the 
outer boundary are included, but are generally small due 
to the high optical depth, and the coarseness of the spatial 
grid considered. 

4. Results 

4--1- Basic molecules: H 2 , CO, and H 2 

Due to their stability, CO and H2O are significant chem- 
ical sinks, with abundances that are relatively constant 
with time. To see this, in Fig. we plot the fractional 
abundance of CO and H 2 throughout the envelope as 
functions of time. As can be seen, the CO abundance is 
essentially constant in time. The abundance has been cho- 
sen to be consistent with observations. 

The water abundance in the warm interior is nearly 
constant, due to the fact that the majority of the oxygen 



Charnley 1997) have run which show that even when the 
oxygen is not initially bound in water, nearly all of the 
available oxygen is converted into water on a timescale 
of about one hundred years due to fast neutral-neutral 
reactions in the warm gas. 

The near discontinuity in the water abundance at 
T ~ 100 K is due to the release of water from grain man- 
tles. This discontinuity is consistent with observations of 
warm (T ~ 300—500 K) water in absorption toward AFGL 
2591 (Helmich et al. |1996| , Boonman et al. [2000D , with the 
lack of strong emission by cold water at long wavelengths 
(Boonman et al. 2000] ), and by detailed modeling of the 
line emisson to be discussed in a forthcoming paper (Boon- 
man et al. - in preparation). 

As noted by Charnley ( |1997 ) , the ion fraction and elec- 
tron density grow with time. As seen in Fig. |^ this leads 
to a destruction of water on timescales of > 10 5 years 
in the interior, in agreement with the results of Charnley 



(1997). While the cosmic-ray ionization continually cre- 
ates ions which destroy water, reformation is temperature 
dependent. A simple extrapolation of the "critical tem- 



perature" for water formation from Charnley (1997) for 
our adopted cosmic-ray ionization rate and density yields 
180 — 200 K. Based upon the temperature structure in Fig. 
|l], this implies destruction of water for r ~ 6-8x 10 15 cm, 
in agreement with the results in Fig. 0. It should be noted 
that the destruction of water for t > 10 5 yrs is probably 
unimportant for AFGL 2591 based both upon the water 
distribution inferred by Boonman et al. (in preparation), 
and the chemical evolution timescale of < 10 5 yrs dis- 
cussed in Sect. 5. 

Finally, the growth in the water abundance with 
time in the exterior occurs through slower (due to the 
lower abundances) ion-molecule reactions. Again, the ion- 
molecule reactions are driven by cosmic-ray ionization. 
In the exterior, average abundances of < 3 x 10 -7 are 
achieved for t ~ 3 x 10 4 years. 

The results in Fig. ^ have interesting implications 
for the interpretation of water abundances. First, a 
simple estimate of the water abundance inferred from 
the model radial column densities [assuming x(H20) = 
N (H2O) / N (H2)] would suggest a fractional abundance of 
water in our model of x(H20) ~3x 10~ 5 . This is a factor 
of 5 lower than the actual water abundance adopted in 
the interior, and would by itself imply a significantly dif- 
ferent structure and chemistry involved. This underscores 
the potential pitfalls in interpreting column densities, as 
well as the importance of modeling the complete physical, 
thermal, and chemical structure of the envelope in order to 
properly compare the relevant regions with observations. 

A second implication is that beam dilution can have an 
important effect on the inferred column densities. A simu- 
lated beam-averaged column density commensurate with 
the beam of the Submillimeter Wave Astronomy Satel- 
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Fig. 3. The radial column densities of H2, CO, and H2O 
as function of time (solid lines). The shaded regions cor- 
respond to the observed abundances (with factor of two 
errorbars). 



Fig. 4. The fractional abundances of HCN, CH4, and 
C2H2 as functions of time for various temperatures. Here 
ft (Ha) = 10 7 cm -3 . Notice the general enhancement of the 
abundances with increasing temperature. 



lite (SWAS) would imply a water abundance of x(H20) = 
N(E 2 0)/N(H 2 ) = 10~ 7 - 10~ 8 depending upon the time 
considered. This low abundance is due to significant beam- 
dilution from the small region of enhanced H2O in the 
large beam. The range of abundances is similar to that 
inf erred by SWAS (se e e.g. Snell et al. |2000| ; Melnick et 



al. 2000 Neufeld et al. 2000 ). Clearly, such an observation 
alone does not constrain the entire envelope. While it im- 
plies that a portion of the envelope (e.g., T < 100 K) has 
a low water abundance, it does not restrict the potential 
for a compact region of significant water abundance. 

As a comparison of the column densities with obser- 
vations, in Fig. H we plot the H2, CO, and H2O column 
densities as a function of time. We assign errorbars of a 
factor of two consistent with the intrinsic uncertainties in 
the H2O and CO results, and with the fact that the H2 
results are scaled from the CO data, as well as various ra- 
diative transfer effects. The ranges of the observed column 
densities are given by the shaded regions. As expected, our 
data match the observed column densities within the un- 
certainties. 

4-2. Hydrocarbon and nitrogen chemistry 



Observations by Lahuis & van Dishoeck (2000) suggest 



that the 14/xm bands of C2H2 and HCN are good trac- 
ers of hot gas. Perhaps more importantly, the increase in 
observed column densities for temperatures above a few 
hundred K implies that their chemistry may be altered 



at high temperatures. Since all of their inferred excitation 
temperatures are well above the expected desorption tem- 
perature of ~ 100 K, it is expected that these enhanced 
abundances are due to warm gas-phase chemistry. 

In order to test this, we have constructed single- 
position models for the chemistry at n(H2) = 10 7 cm~ 3 
and T > 200 K. The results are shown in Fig. [|. Clearly, 
higher temperatures do increase the abundances of simple 
hydrocarbons and nitrogen-bearing species, with higher 
abundances prevalent once T ~ few hundred K. 

The enhanced HCN abundance is similar to that found 
by Rodgers & Charnley ( 200l| ). In parallel with their work, 
we find that the 756 K endothermic reaction CN + H2 — > 
HCN proceeds quickly for T > 200 K, producing signif- 
icant HCN. However, while Rodgers & Charnley fl2001| ) 
assume the reaction C++NH 3 favors HCNH + (following 
ab initio calculations by Talbi & Herbst 1998), we assume 



that H2NC + is the favored product to account for the ob- 
served HNC/HCN abundance ratio in many sources. In 
our case, then, the CN is formed via the neutral-neutral 
reaction CS + N — > CN + S. This reaction has a barrier of 
1160 K, leading to significant production for temperatures 
above 200 K. Overcoming these barriers can increase the 
abundance from a peak of 10 -8 at 200 K, to ~ 10~ 7 for 
t > 4 x 10 4 years, and ~ 10~ 6 for t > 3 x 10 5 years for 
T > 400K. 

Methane shows perhaps the most dramatic increase in 
abundance at very high temperatures. In fact, methane 
contains more carbon at the latest times than all species 
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other than CO2 and CO initially. This is due to the fact 
that ion-molcculc reactions driven by cosmic-ray ioniza- 
tion (e.g., He+ + CO -» C+ + O) can pr oduce C+. This 
then reacts via carbon insertion (Herbst 1995) with H2 
to form CH + at high temperatures, and then in a chain 
with H2 up to CHj, which dissociatively recombines to 
form CH. While CH^ can also dissociatively recombine to 
form CH2, the dominant pathway to CH2 at high temper- 
atures is CH + H 2 + 1760K -> CH 2 + H. Reactions with 
H 2 then produce CH4 (overcoming barriers of 6400 K and 
4740 K to form CH3 and CH4 respectively), leading to 
abundances of 1 — 3 x 10~ 7 for T < 400K. However, once 
the temperature increases to ~ 600 — 800K, abundances 
can reach x(CH4) ~ 10~ 6 at t ~ 3 x 10 4 years. 

Acetylene is also enhanced at high temperatures. The 
pathway here is similar to that in diffuse and dark clouds 
(van Dishoeck & Hogerheijde 1999). However, in our 
model, acetylene is formed via reactions of water with 
C 2 H^" instead of dissociative recombination. A second dif- 
ference is that C 2 H 3 f " is produced via C 2 H 4 + H+. While 
the "usual" CH4 + C + production route still occurs, the 
destruction of C 2 H4 by O is reduced as the temperature 
increases due to the fact that the oxygen is quickly con- 
verted into water by neutral- neutral reactions (see Sect. 
4.1). Again, cosmic-ray ionization, carbon insertion, and 
water play a role, both in the production of H;|~, and in the 



production of C 2 H4 via CO — * C^ 



CH++CH4 



C 2 Hg + H 2 — > C 2 H4. The enhanced C 2 H2 abundance 
is in the range 5 x 10~ 9 < .x(C 2 H 2 ) < 2 x 10~ 8 for 
200K < T < 800K at t = 3 x 10 4 years. At late times, 
it is almost always less than 3 x 10~ 8 at 200 K, less than 
5 x 10~ 8 at 400 K, and can reach 5 x 10~ 7 at 800 K. 

In order to see how this high-temperature chemistry 
pertains to our model, in Fig. || we plot the fractional 
abundances of HCN, CH4, and C2H 2 throughout the enve- 
lope for various times. As expected from the previous dis- 
cussion, we see enhanced abundances of HCN, C2H 2 , and 
CH4, especially in the warm interior. The enhancement 
of C 2 H 2 in the exterior has two primary causes. First, in 
this region C 2 H 2 is primarily formed via dissociative re- 
combination of C 2 H^". The destruction of C 2 H^" by O has 
a 215 K barrier that cannot be overcome in the cool ex- 
terior, leaving more C 2 H^" to produce acetylene. Second, 
an alternate production pathway via C3H3 +0— >C2Ha is 
enhanced in the exterior due to our increased initial O 
abundance in that region (see Table |). 

Cosmic-ray driven ion-molecule chemistry again plays 
a role for t > 10 5 years. In particular, the destruction 
of HCN near 10 16 cm is due to reactions with HCO + . 
For C2H 2 both HCO + and O are important destruction 
reactants near 10 16 cm. The enhancement in C 2 H 2 near 
r~6-8x 10 15 cm is due to a decrease in atomic oxygen 
at this position for late times (see also Sect. 4.6). 

Observations of high- lying HCN lines in the submil- 
limeter were undertaken by Boonman et al. ( 2001 ). They 
utilized a sophisticated radiative transfer model of the ex- 
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Fig. 5. The time evolution of the fractional abundances of 
HCN, CH4, and C2H2 throughout our model, incorporat- 
ing the temperature and density distributions desribed in 
the text. The HCN data are labeled by the time in years, 
where a(b) = a x 10 b . The times for each curve increase 
upward at the inner radial position. Note the enhancement 
of the abundances in the warmer interior. 



citation, line shapes and strengths to analyze their data, 
and suggested that HCN follow a "jump" model, with an 
abundance of x(HCN) - 1 x 10~ 6 for T £ 230K, and 
x(HCN) ~ 10 -8 for the cool exterior. The results in Figs. 
|] and [5] are consistent with this supposition, with abun- 
dances of a few x 10~ 7 at high temperatures, and ~ 10 -8 
at lower temperatures and later times. 

As expected, the results in Fig. || are not as dramatic 
as in Fig. Q as our physical and thermal model only ex- 
tends into T ~ 440 K, less than the temperatures at which 
the greatest enhancements occur. Consequently, care must 
be used when comparing the results with observations, as 
the different temperature components may not necessarily 
probe the portions of the region being modeled. 

Such a comparison is given in Fig. |^. Here the model 
predictions for CH 4 , C2H2, and HCN are compared with 
the infrared observational data, which probe column den- 
sity. In the left-h and panels we compare to the data of 
Carr et al. ( 1995 ), omitting the 200 K data as these arise 
in the outflow. In the right-hand panels we compare to the 
data of Boonm an (p rivate communication) , and Lahuis & 
van Dishoeck ( |2000| ). 

When we compare with the lower temperature data, 
the CH4 model results are close to the observed error 
bounds. On the other hand, they are well above the high 
temperature component of the observations. This is not 
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component fit by Lahuis & van Dishoeck (2000). In both 



Fig. 6. A comparison of the predicted and observed col- 
umn densities of HCN, CH4, and C2H2 as a function of 
time. The model predictions are given by the solid lines 
and accompanied by the filled circles. The observations are 
divided into two groups. The left-hand panels are for the 
infrared data of Carr et al. ( 1995 CELZ), while the right- 
hand panels are for ISO data from Lahuis & van Dishoeck 



(2000 Lv) and Boonman (B(pc), private communication). 
Data which are upper limits are signified by downward 
arrows. Other data have been given an arbitrary factor of 
3 uncertainty, and are given by the shaded regions. 



a suprise, as the CH4 chemistry is relatively unaffected 
below about 400 K, with very significant production at 
higher temperatures. 

On the other hand, the C2H2 data fit the high tem- 
perature components of the observations. This seems to 
imply that while high temperature chemistry can be im- 
portant, the effects are noticably smaller than for CH4, 
consistent with the results of Fig. [|. In particular, the 
fact that the predicted column density is so much higher 
than the low-temperature column density suggests that 
warm chemistry can enhance C2H2, while the fact that 
the predicted column densities fall in the lower range of 
the observed values suggests that there exists room for 
some enhancement (~ 3 — 5x) in the C2H2 abundance at 
higher temperatures, consistent with Fig. |5[ 

Finally, in the lower panels of Fig. ^ we show the com- 
parison for HCN. Here we see the potential for further im- 
portance of high-temperature chemistry. In the lower-left 
panel, the HCN model prediction is consistent with the 
upper limit derived by Carr et al. ( 1995| ) at low temper- 
atures. Similarly, in the lower-right panel, the predicted 
column densities are consistent with the low-temperature 



cases, it appears that our model reproduces the produc- 
tion of cool HCN quite well. 

On the other hand, the model predictions are well be- 
low the observed column densities for the hot components 
in each of the panels. This is most probably due to the 
significant production of HCN at temperatures above 400 
K (see above). This is further supported by the fact that 
when a single temperature component is determined by 
Lahuis & van Dishoeck ( |2000[ ), they find T ~ 600 K. Taken 
at face value, their data suggest that our model does not 
extend inward far enough to include this hot gas. 

At this point, one may ask if a simple extension of our 
power-law model inward would increase the temperature 
and column density sufficiently to fit the observed HCN 
data (i.e., at 1010K). We have examined this possibility 
by extending our model inward, with no success. While a 
fractional abundance of x(HCN) ~ 10~ 7 would reproduce 
the data, the conditions necessary would also produce a 
water column density A^(H20) ~4x 10 19 cm~ 2 , over an 
order of magnitude above the observations. 

An alternative solution is to adopt a "flattened" (i.e., 
n(r) oc r°) density profile for r < r in . In this case, the 
extra column of water would be consistent with the obser- 
vations, and the column of HCN would vary as iV(HCN) ~ 
2 x 10 15 (x(HCN)/10- 7 ). While the column could be fit if 
x(HCN) = 10~ 6 , this is inconsistent with the results of 
Fig. ^. First, the chemistry does not show strong varia- 
tion between 400 and 800 K, suggesting that high tem- 
peratures alone will not produce significantly more HCN. 
Furthermore, to achieve x(HCN) = 10~ 6 at these tem- 
peratures would require an extended time for chemical 
evolution in the interior, and would be inconsistent with 
the abundances of other observed species (see Sect. 5). 

There are four possible resolutions to this difficulty. 
First, and least likely, is the possiblity that the chemical 
evolution time in the interior is somehow longer than in 
the exterior. We can think of no way in which this may 
occur. The second possibility is that the hydrocarbon and 
nitrogen chemistry is currently incomplete, especially at 
high temperatures. If another pathway to producing HCN 
exists above about 600 K, it would be possible to have 
abundances of 10~ 6 . Third, there is the possibility that 
HCN is present in grain mantles, and is injected into the 
hot gas. Thou gh th is is expected to be unimportant (van 
der Tak et al. 1999 ), it may conceivably play a small role. 

The fourth, and perhaps most likely, possibility is that 
there exists some as of yet unidentified destruction mech- 
anism for water at high temperatures. This would remove 
the problem of the overly-large water column if the en- 
velope were to simply extend further inward. It is pos- 
sible that evidence exists for this. As discussed by van 
Dishoeck ( |1998| ) observations of water gas and ice toward 
various sources show significantly less total water in hotter 
sources than in cooler sources. Given our current under- 
standing of the chemistry of H2O production, it would be 
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Fig. 7. The fractional abundance of CS throughout the 
envelope for various times. Nearly half of the sulphur is 
in CS at late times in the cool exterior, essentially "fix- 
ing" the gas-phase sulphur abundance. The agreement 
with observations in the warm interior, however is not 
fixed. The curves are labeled by the time in years, where 
a(b) = a x 10 h . 



Fig. 8. The fractional abundance of SO2 throughout the 
envelope at various times. Note the increase near T — 
100 K as the free sulphur is forced into SO2. The decrease 
at high temperatures is due to the loss of the reactant 
OH via its more efficient inclusion into water at those 
temperatures. The curves are labeled by the time in years, 
where a(b) = a x 10 b . 



easiest to explain this effect if there existed a mechanism 
for H2O destruction at high temperatures. Further study 
into the high temperature chemistry of water, hydrocar- 
bons, and nitrogen-bearing species would be of significant 
importance in understanding this problem. 

4-3. Sulphur chemistry 

The chemistry of sulphur in hot cores is well-described by 
Charnley p97)) . In our model, we have adopted a chem- 
istry and set of initial conditions (in the warm region) 
which is similar. However, given the fact that his model 
was for a single point in space, while our model extends 
over a range of physical and thermal parameters, and given 
recent observations of sulphur-bearing molecules toward 
AFGL 2591, we present our results here. 

In the cool exterior of our model we find that there 
exist a large number of pathways to shuttle sulphur into 
CS. The end product is that approximately 50% of the 
sulphur is transformed into CS by t ~ 10 5 years. This is 
shown in Fig. ^, where we plot the fractional abundance 
of CS for various times. No single production reaction ac- 
counts for more than 25% of the final CS abundance. This 
means that, at late times at least, CS is a good measure of 
the sulphur abundance in the exterior. To accomodate this 
fact, and in order to match observations of the CS abun- 



dance (see Table [j]), we adjust the initial sulphur abun- 
dance to x(S) = 6 x 10" 9 for T < 100 K. This produces a 
nearly constant abundance in the exterior in good agree- 
ment with the observations. 

In the interior, the CS abundance increases at inter- 
mediate and late times to x(CS) ~ 10~ 8 . This is also 
in agreement with the observations. However, while the 
abundance in the exterior is essentially "forced" by our 
initial sulphur abundance, the fraction in CS in the inte- 
rior is not. 

The variation in the interior CS abundance (as with 
water and the hydrocarbons) is again related to the oxygen 
and cosmic-ray driven ion-molecule chemistry. In particu- 
lar the dip near r~5x 10 16 cm is due to the increased 
atomic oxygen abundance (see Sect. 4.6) in this region. 
This leads to more conversion of sulphur to SO2, and thus 
less to CS. The decrease in CS abundance near 10 16 cm 
at t = 1 — 3 x 10 5 yrs is due to the fact that there is less 
OH available for conversion of sulphur out of H2S to CS. 

The sulphur abundance in the warm (T ~ 100 — 400 K) 
gas is well-determined by the SO2 abundance. In our 
model, S0 2 is formed by H 2 S + (OH,H) -> HS + O -> 
SO + OH -> S0 2 . The initial reactions of H 2 S with H 
and OH have barriers of 352 K and 80 K, respectively. 
As a result, little SO2 is produced in the cool exterior, 
while the barriers can be overcome in the interior leading 
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to significant SO2 production. As the temperature further 
increases, however, the OH can be more easily forced into 
water, leaving little for the SO + OH — > SO2 + H reaction. 
This can be seen in the very interior of Fig. |8[ where the 
SO2 abundance drops at high temperatures. In our model 
approximately 90% of the sulphur returns to atomic form 
at ~ 440 K, with approximately 10% in H2CS, and a few 
percent in CS and OCS. 

While we are are unable to identify the sulphur resevoir 
assuming solar abundances roughly hold, it appears that 
a significant portion would need to exist in or on dust 
grains. Under this constraint, we can also identify SO2 
as the primary sink of molecular sulphur in warm (100- 
300 K) gas (assuming no O2 is released during heating of 
the grain mantles - Charnlcy 1997). As a result, the sul- 
phur abundance in warm molecular gas at later times can 
be approximately determined by the SO2 abundance. In 
our model, this requires the adjustment of the initial H 2 S 
abundance from the value of 10 -7 adopted by Charnley 
( |1997| ) to 1.6 x 10~ 6 . This value is, coincidentally, s imilar 
to the H2S gas-phase abundance seen by Minh et al. ( 1990 ) 
toward Orion. A comparison o f our m odel predictions with 
observations by Keane et al. ( 2001 ) show similar column 
densities of 4 x 10 16 cm -2 and 6 x 10 16 cm -2 respectively. 
It is also intruiging that the excitation temperature in- 
ferred by Keane et al. ( |200l| ) for S0 2 toward AFGL 2591 
is ~ 750K, suggesting formation in a warm dense region 
of a few hundred K. 
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Fig. 9. The column density of CO2 predicted by our base 
model as a function of time (solid and dashed lines) for two 
different assumed desorption temperatures. No impulsive 
heating event is assumed (see text). The observed column 
is shown by the shaded region. Notice the extent to which 
the model overpredicts CO2. 



4-4- CO2 chemistry: potential heating events? 

An important problem in the chemistry of the envelopes of 
massive young stars is the low observed gas-phase abun- 
dance of CO2 (see e.g., van Dishoeck & van derTak ^OOO]) . 
Observations by I SO in dicate large soli d CO 2 abundances 
(de Graau w et al. |1996| , W hittet et al. |l998j Ehrenfreund 
et al. |1998| , Gerakines et al. |l999| ), with a C0 2 /H 2 abun- 
dance in the ice mantles of 10 — 20%. In the warm regions 
close to the protostars, these mantles should be evapo- 
rated. Assumi ng w ater ice abundanc es of a few xlO -5 
(Tielens et al. 1991; Gensheimer et al. 1996) implies a lib- 
erated fractional abundance of x(G02) ~ 10 -5 — 10~ 6 . 
On the other hand, ISO observations of gas-phase CO2 
(van Dishoeck et al. 1996, Boonman et al. 200C) suggest 
x(G02) ~ 10~ 7 . These results indicate that CO2 is quickly 
destroyed after evaporation from ice mantles. 

To see this discrepancy between the amount of CO2 
predicted in our base model and that observed, in Fig. || 
we plot the predicted CO2 column density as a function of 
t ime. Also plotted are the observations of Boonman et al. 
( 200C ). Clearly, the base model significantly overpredicts 
the CO2 column density in AFGL 2591, confirming the 
general results above. 



Charnley & Kaufman (2000) studied destruction of 
CO2 by both H and H2, suggesting that destruction of 
CO2 by H in postshock flows could be important. In or- 



der to test this, we have constructed models with non- 
zero atomic hydrogen abundances, as would be expected 
in partially dissociative shocks. While the CO2 can be ef- 
fectively destroyed on a shock cooling timescale of ~ 30 
yrs, CH4, NH3, and H2O can be destroyed even more effi- 
ciently. While this does not pose a significant problem for 
CO2 or NH 3 which have low observed abundances or upper 
limits, there are effects on other species. In particular, in 
dissociative shocks the water column density is decreased 
by a factor of 2-3. Furthermore, once destroyed, only lit- 
tle water is re-formed in the range 100 < T(K) < 300, 
inconsistent with the results of Boonman et al. (in prepa- 
ration). Similarly, the O and O2 abundances are signifi- 
cantly increased. On the other hand, the CH 4 abundance 
is decreased by an order of magnitude in the interior. This 
process only requires a few percent H 2 dissociation. 

A second potential difficulty is that it is unclear if a 
large enough fraction of the envelope can be disturbed by 
a shock to significantly affect the global CO2 abundance, 
as evidenced by the relatively small line-widths in much 

1999|) . 



of the envelope (van der Tak et al. 



Doty et al. (2002) reconsidered this problem in light 



of previously unused laboratory measurements of the de- 



struction of CO2 by H 2 (Graven & Long 1954). They found 



that destruction by H2 may dominate destruction by H in 
the very warm gas, near T ~ 1000 — 1600 K. Whil e this 
may occur in a number of ways, Doty et al. ( |2002| ) con- 
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sidered two possibilities: a uniform temperature increase 
(such as from the passage of a v ~ 20 — 30 km s _1 MHD 
shock - Draine et. al. 1983), and a central luminosity in- 
crease caused, for example, by an accretion (FU-Orionis- 
type) event. The possibility of impulsive heating events 
may be supported by evidence from continuum emission 
by crystalline silictes (Smith et al. 2000, & Aitken et al. 
1988| ) which suggests that an annealing event may have 



occurred in AFGL 2591. If such a heating event occurred, 
Doty et al. ( 2002 ) find that it is possible for the CO2 to 
be removed on a timescale of 10° — 10 4 years by H2. Re- 
cent calculations of the potential surface for the CO2 + 
H2 reaction suggest, however, that the barrier for the re- 
action may be higher than indicated by the old laboratory 
experi ments, so that this issue remains unsettled (Talbi & 
Hcrbst 2002 ). Clearly, further laboratory studies of this 
reaction at high temperatures are urgently needed. 

While speculative, destruction of CO2 by H2 in this 
fashion has some advantages. First, there is very little 
atomic hydrogen available to affect the chemistry, and in 
particular to influence CH 4 , 0,02, and H 2 0. Second, and 
perhaps more importantly, variations in the observed col- 
umn density of CO2 may potentially be explained by vari- 
ations in the size and/or duration of the proposed heating 
event - depending upon its origin, or the time since the 
heating event and the local cosmic-ray ionization rate. 

As a final note, it is interesting to also consider the 
possibility that the CO2 desorption temperature may be 
greater than 100 K. Recent work by Fraser et al. (2001), 



suggests that the desorption temperature of water may be 
as high as 120 — 130 K. If the solid CO2 is contained in a 
water- ice m atrix as suggested by observations (Gerakines 
et al. 199E ), then it may be interesting to consider the 



effect of this higher desorption temperature on iV(C02). 
In Fig. |], we present predicted column densities for the 
re-formation of CO2, assuming desorption temperatures 
of both 100 K and 130 K. The effect is a decrease in the 
CO2 column densities by a factor of two, insufficient to 
explain the discrepancies. 

4-5. Cosmic-ray ionization rate 

As discussed earlier, cosmic-ray ionization can play an 
important role in driving ion-molecule chemistry at later 
times. In our model, we adopt the cosmic-ray ionization 
rate for AFGL 2591 of C = 5.6 x 10 



-17 



1 as determined 
While the cos- 



by van der Tak & van Dishoeck (|2000j) 
mic ray flux is unique, the ionization rate will vary with 
position if the particles are absorbed. As evidence for cos- 
mic ray absorption is inconclusive (see, e.g., van der Tak 
200 1 ), we adopt a single cosmic ray ionization rate for 
AFGL 2591. 

In Fig. we plot the predicted fractional abundance 
of HCO + and N2H+. There are two important features. 
First, there is significant destruction of HCO + at the water 
desorption position, due to the reaction HCO + + H2O — > 
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Fig. 10. The fractional abundance of HCO + and N 2 H+ 
througout the envelope for various times. Note the marked 
decrease about 100 K, where reactions with H2O become 
important. The curves are labeled by the time in years, 
where a(b) = a x 10 b . 

H30 + + CO. This is in agreement with the model of van 
der Tak & van Dishoeck (|2000] ). While they argue that this 
jump in abundances is not important in constraining the 
cosmic ray ionization rate, our overall HCO + abundance is 
consistent with their observations, and thus lends support 
to their somewhat high value for £ in AFGL 2591. The 
situation is similar for N2H" 1 ". Second, at t = 3 x 10 5 yrs, 
the ion abundances increase in the interior. This is con- 



sistent with Charnley (1997), and is due to the fact that 
the cosmic-ray ionization continues to produce more ions, 
which eventually destroy a significant fraction of the com- 
plex molecules up to the position where the temperature 
is high enough to re-form them. 

The cosmic-ray ionization rate also affects the abun- 
dance of ELt. In our models, reasonable time (i > 10 3 
years) column densities are ~ 2.6 x 10 13 cm~ 2 , almost an 
order of magnitude below those observed by McCall et al. 
( gggg ). If a comparison of these results were used to in- 
fer a cosmic-ray ionization rate, one would obtain a much 
larger value. While large H^" abundances in the diffuse 
ISM have been reported by McCall et al. ( |2002j ) - which 
they suggest may be due to uncertainties in disso ciative 
recombination rate - van der Tak & van Dishoeck (2000) 



have also noted there exists a variation in Hjj~ column den- 
sity with distance which suggests that intervening clouds 
may be important. 

The cosmic-ray ionization rate also affects the HCN 
abundance. Decreasing £ makes it harder to form HCN. 
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Fig. 11. The fractional abundance of 0, OH, and O2 
throughout the envelope for various times. The curves are 
labeled by the time in years, where a(b) = a x 10 b . The 
times increase upward. The only exception is the dip for 
OH near 10 16 cm, which corresponds to t = 3 x 10 5 yrs. 



For example lowering £ by a factor of three decreases the 
enhancement of HCN by a factor of three even at 800 K, 
placing the warm HCN abundance at < 3 x 10~ 7 - well be- 
low the observations. Furthermore, the same change also 
increases the time for the cold column of HCN to reach 
the observed range to t ~ 10 5 yrs, in disagreement with 
the age constraints discussed below in Sect. 5. 

Based upon these results, it appears that the value 
of £ = 5.6 x 10~ 17 s _1 inferred by van der Tak & van 
Dishoeck ( [2000 ) is correct to within a factor of three. Any 
value much lower would significantly hamper the produc- 
tion of HCN, making for disagreement with the observa- 
tions. Any value much higher would be in conflict with 
the observed ion abundances. 



4-6. Other species 

As oxygen and oxygen-bearing species can have a signif- 
cant effect on the chemistry, in Fig. [II] we plot the frac- 
tional abundances of O, OH, and O2 as functions of posi- 
tion for various times. 

The increase in the atomic oxygen abundance near 
10 16 cm is due to the fact that O is freed from water at 
late times via ion- molecule reactions with H2O and CO as 
discussed in Sect. 4.1. As the water is destroyed, the main 
production mechanism for OH [(HCO + , H3 ) + H 2 — > 
H30 + +c~ — » OH] is removed. This leads to a decreased 
OH abundance at this position at late times. 



The peak in the atomic oxygen abundance near r ~ 
5 — 8 x 10 15 cm is due to the competition between pro- 
duction of O by ion- molecule reactions with CO, and the 
destruction of O at high temperatures by reactions with 
OH and H2. Once the temperature reaches ~ 180 — 200 K, 
neutral-neutral re-formation of water can balance the de- 
struction by ion-molecule reactions on these timescales, as 
discussed in Sect. 4.1, and in Fig. ^J. This leads to a greater 
OH abundance at these positions, and thus a decreased O 
abundance. 

In any case, the excess atomic oxygen is easily con- 
verted to molecular oxygen over time at temperatures less 
than 300 K. This places an important constraint on the 
temporal evolution of the source as discussed in Sect. 5 
below. 

It is also interesting to note that the dominant nitrogen 
resevoir is molecular nitrogen. While atomic nitrogen is 
somewhat abundant (see Table ||), only about 1 % or less 
of the nitrogen is in atomic form - and that preferentially 
at later times. 

Although we have endeavored to consider detailed 
comparisons between our model predictions and obser- 
vations, a worthwhile test of any model is the predic- 
tions it makes for future observations. Consequently, in 
Table S, we give predicted radial and beam-averaged col- 
umn densities at t = 3 x 10 4 yrs for various species with 
N > 10 13 cm -2 . The beam-averaged column densities as- 
sume a gaussian beam of full-width at half-max of 15 arsec, 
though the results are insensitive to this assumption. 

4-7. Implications for mantles and mantle destruction 

Finally, although we do not explicitly consider grain- 
surface chemistry, it is worthwhile to discuss the impli- 
cations our results have on the grain mantles, and grain- 
surface chemistry. Two observed species that may form on 
grain mantles are H2CO and CH3OH. In Fig. [l2] we plot 
the fractional abundance of H2CO and CH3OH through- 
out the envelope at various times. 

When low-temperature depletion is assumed, the col- 
umn densities for H2CO and CH3OH are generally close 
to the observations. In particular, while van der Tak, 



van Dishoeck & Caselli (200C) report column densities of 
N(R 2 CO) = 8 x 10 13 cm^and W(CH 3 OH) = 1.2 x 10 15 
cm -2 respectively, the model predicts a CH3OH column 
density about 6 times lower, and an H2CO column density 
about 5 times higher. 

On the other hand, detailed radiative transfer model- 
ing by van der Tak, van Dishoeck, & Caselli ( |2000 ) sug- 
gests that the observed lines are consistent with a uniform 
H2CO abundance of 4 x 10 -9 , and a CH3OH abundance 
of 2.6 x 10~ 9 for T < 100 K, and 8 x 10~ 8 for T > 100 K. 
For comparison the predicted abundances for these species 
are shown in Fig. O. The H2CO abundance in the cool 
exterior is consistent with the inferred abundance, while 
the abundance in the warm interior is predicted to be 
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Table 3. Predicted Column Densities at t — 3 x 10 4 yrs. 



Species 


Radial (X) 


iVbcam(X) - 15" 


OI 


4(18) 


3(18) 


N 


3(16) 


2(16) 


S 


1(16) 


7(13) 


NO 


1(16) 


5(15) 


OH 


9(15) 


1(15) 


SO 


4(15) 


2(13) 


H 2 CS 


2(15) 


1(13) 


C 3 H 


2(15) 


2(15) 


C 4 H 


1(15) 


1(15) 


C3H2 


7(14) 


7(14) 


CH3OCH3 


6(14) 


7(12) 


CHOOH 


4(14) 


1(14) 


NH 2 


3(14) 


2(13) 


CH2CO 


2(14) 


2(14) 


CH 3 


2(14) 


6(12) 


NH 


1(14) 


3(12) 


H3O+ 


8(13) 


2(13) 


CN 


5(13) 


5(13) 


OCN 


5(13) 


1(13) 


NS 


5(13) 


9(10) 


C 2 S 


4(13) 


3(13) 


HS 2 


4(13) 


1(12) 


C 6 H 


4(13) 


4(13) 


C3H3 


3(13) 


2(13) 


CCN 


3(13) 


3(13) 


H2C3 


3(13) 


2(13) 


CH3OH+ 


2(13) 


2(11) 


HNO 


2(13) 


1(13) 


CH3CHO 


2(13) 


7(12) 



a(b) means 
All column 



a x 10° 

densities given in cm" 2 



significantly higher and decreases only slowly with time. 
The CH3OH, on the other hand, does not fit the inferred 
abundances very well. While there does exist a "ju mp" a s 
suggested by van der Tak, van Dishoeck, & Caselli fl2000] ), 
the abundances predicted by the model are significantly 
too low in the exterior and too high in the interior. 

For comparison, we have also run models where the 
abundances of H2CO and CH3OH are initially undepleted 
for T < 100 K. In models where the cold initial abun- 
dances of CH3OH are set equal to the hot initial abun- 
dances only CH3OH, C3H3, and CH3OCH3 show column 
density differences of a factor of three or more. When 
we adopt this initial abundance, the column density of 
CH3OH increases to 4.4 x 10 15 cm" 2 , a value approxi- 
mately 4 times larger than the observations. These re- 
sults suggest that while we can reproduce the CH3OH col- 
umn density with some accuracy, simple gas-phase chem- 
istry alone cannot reproduce the apparent details of the 
CH3OH abundance distribution. As a result, it appears 
that CH3OH can be strongly affected by grain surface 
chemistry, also in the cooler regions. 
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Fig. 12. The fractional abundance of H 2 CO and CH 3 OH 
throughout the envelope for various times. The curves are 
labeled by the time in years, where a(b) = a x 10 b . Note 
that both species are initially depleted from the gas phase 
for T < 100 K. 
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Fig. 13. The fractional abundance of H 2 CO and CH3OH 
throughout the envelope for various times, after a hypo- 
thetical non fiVdissociative heating event. The curves are 
labeled by the time in years, where a(b) = a x 10 b . Note 
that both species are initially depleted from the gas phase 
for T < 100 K. 
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Allowing a cold H2CO abundance equal to the warm 
abundance changes the predicted column density by only 
a factor of two, with only minor differences for all other 
species. This implies that only observations of high-enough 
spatial resolution to differentiate between the warm and 
cold phases, or use of high-excitation lines will be able to 
best determine the nature of H2CO formation. 

It is also interesting to consider the fact that submil- 
limeter observations suggest that the abundances of warm 
H 2 CO and CH3OH are factors of 100 - 1000 below the 
solid state abundances. In our models, we find that the 
ratio of warm to cold H2CO for a 15 arcsecond beam 
is 60 < [Af(H2CO)]T>ioo/[A r (H 2 CO)]T<ioo < 500 for 

3 x 10 4 < t(yrs) < 10 5 . On the other hand, for CH 3 OH we 
find that 1 < [A^(CH3OH)] T > 100 /[A^(CH 3 OH)] T < 100 < 2 
for 3 x 10 4 < t(yrs) < 10 5 . This confirms the previ- 
ous suggestions that that gas-phase chemistry may per- 
haps dominate grain-surface chemistry in the production 
of H2CO, while there must be some other (presumeably 
grain-surface) pathway to the production of CH 3 OH. 

As a final note, we have considered the effect of a heat- 
ing event (as proposed in Sect. 4.4) in which H2 is not 
dissociated on the H 2 CO and CH 3 OH chemistry. The re- 
sults are shown in Fig. 13. As can be seen, the later-time 
abundances are more consisent with the results inferred 
by van der Tak, van Dishoeck, & Caselli (|200C| ). In partic- 
ular, the CH3OH abundance in the interior is in the range 
of 5 x 10~ 8 < a;(CH 3 OH) < 6 x 10~ 7 , while in the exterior 
the abundances can reach as high as 0.3 — 1 x 10~ 9 . Like- 
wise, the H 2 CO abundance for 3 x 10 4 < t(yrs) < 3 x 10 5 
is in the range 10" 7 > x(H 2 CO) > 3 x 10~ 10 . While not 
conclusive, this brackets the inferred H2CO abundance of 

4 x 10~ 9 nicely. If further suggestions of a heating event 
are found, it may be useful to re-visit these data for com- 
parison with observations as they may provide a gas-phase 
mechanism for the production of H2CO and CH3OH. 



5. Time constraints 

Based upon the large amount of observational data for 
AFGL 2591 (see Table |Tj), and given the time-dependent 
nature of the reaction network, one important test of the 
physical and chemical model would be a determination 
of the chemical evolution time of the envelope, consistent 
with all or most of the observed species. This has been 
proposed and carried out previously (e.g., Stahler 1984; 



Millar 199C; Helmich et al. 1994; Hatchell et al. 1998) in 



single-point models of dense cloud cores, with some suc- 
cess. 

In our case, we determine the time-dependent frac- 
tional abundances and column densities for each of the 
species observed in Table [ij As discussed in Sect. 2.3, we 
divide the data into two sets: those for which infrared 
absorption measurements have been made, and those for 
which submillimeter emission measurements have been 
made. 
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Fig. 14. A comparison of the predicted and observed 
abundances and column densities for the species listed in 
Table [l], and observed in the infrared (see text) . The solid 
lines correspond to agreement between the models and 
observations within a factor of 3, and the dashed- lines to 
within a factor of 10. The species are listed, with notes on 
the observational fits given as parentheses as in Table [jj 
The two shaded regions denote the regions of potential and 
preferred fit between the model and the observations (see 
text). Notice the agreement with Fig. [L5L 



In Figs. |l4| and |l5] we plot the approximate time ranges 
over which the listed species match the observed data. In 
both figures, the solid lines represent agreement to within 
a factor of 3, while the dashed lines represent agreement 
to within a factor of 10. These limits can be considered 
good and acceptable le vels o f agreement respectivel y (see 
e.g., Millar & Freeman 1984 and Brown & Charnley 1990 
In both figures, the lower limits to the chemical evolu- 
tion time are capped at 10 3 years. As in Table [j], observa- 
tional data are appended to each species name. We include 
all species from Table [l], except for CO2. As discussed in 
Sect. 4.4, there are significant discrepancies and questions 
about the gas-phase chemistry of CO2, and as such we 
have treated it separately in that section. For species other 
than CO2, we include the data that is relevant to compar- 
ison with our model (i.e., the most reliable components). 

In Fig. [l4|, we can see a wide variation of possible times 
when constraining the models by the infrared data. We 
place two limits on the evolution: a wide range of 10 3 < 
f (yrs) < 1 x 10 5 , and a preferred limit of 3 x 10 3 < t(yrs) < 
2.5 x 10 4 . These regions are identified by the shading in 

Fig. n 
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Fig. 15. A comparison of the predicted and observed 
abundances and column densities for the species listed in 
Table [l], and observed in the submillimeter (see text). The 
solid lines correspond to agreement between the models 
and observations within a factor of 3, and the dashcd-lines 
to within a factor of 10 (see text for details). The species 
are listed, with notes on the observational fits given as 
parentheses as in Table § The shaded region denote the 
regions of preferred fit between the model and the obser- 
vations. Notice the agreement with Fig. [l4|. 

For comparison, in Fig. |lj we plot the constraints on 
the time for the submillimeter data. As discussed in Sect. 
2.3, results based upon sophisticated, self-consistent, ra- 
diative transfer modeling are given slightly more weight, 
and identified by the bold lines in Fig. [TBI Also, the O2 
upper limit by SWAS is given more weight. This is done 
because in the absence of upper limits for the O2 abun- 
dance toward AFGL 2591 in particular, we have used the 
largest quoted up per li mit of N(0 2 )/N(H 2 ) < 9 x 1CT 7 
(Goldsmith et al. 2000] ). Where radiative-transfer model- 
ing derived abundances are not available, we have calcu- 
lated the appropriate beam-averaged column densities for 
comparison with the observations. 

Most species in Fig. [l5] fit the models to within an or- 
der of magnitude of the observational data. Those that do 
not fit are not significant defects, for a number of reasons. 
First, not all data are in disagreement with the models - 
both OCS and HC3N have other observations / reductions 
which do agree with the models. Second, the discrepancies 
can be understood on a case-by-case basis. For instance, 
the chemistry and reaction rates of OCS are only poorly 
understood at best (Millar, private communication). It is 
interesting to note, however, that the radial OCS column 



density matches the observed column density. Also, we ex- 
pect the potential difficulties with species related to HCN 
(such as HC3N and CH3CN) as our model does not probe 
the complete region over which significant HCN produc- 
tion may be important (see Sect. 4.2 above). In the case 
of HCS + , the abundance is strongly affected by enhance- 
ments of CS abundance at temperatures of 10 — 20 K 



(Helmich 1996), lower than all but our outermost tem- 



perature, signifying that the AFGL 2591 envelope may be 
more extended than we have assumed. In a similar fash- 
ion, CN is str ongly influenced by UV radiation from a 
PDR (Helmich |l996 ), a radiation source not considered in 
our model. 

The age-constraints implied by the results in Fig. [l^ 
suggest chemical evolution times in the range 7 x 10 3 < 
t(yrs) < 5 x 10 4 , with a strong preference for t~3x 10 4 
years. These constraints are shown by the shaded region 
in the figure. 

In order to attempt to quantify this result, in Fig. 
|l6| we plot the chi-squared value between the models 
and the observations, normalized to the maximum chi- 
squared in the entire time evolution. The x 2 is defined by 
X 2 = Ejiwf x [j/modou - Vohsd] 2 1 °f ■ Here we include the 
weight from Table |as m, and assume uncertainties of a 
factor of 5 in the observations for all species except O2, 
for which we assume a factor of 2 uncertainty as we al- 
ready have adopted the highest observed upper limit from 



Goldsmith et al. (2000). These results are generally con- 
sistent with those of Figs. [l4| and [l5| both in terms of the 
preferred times as well as the relatively lower level of con- 
straint provided by the IR data. In particular, while times 
of up to ~ 10 yrs may be possible, there appears to be a 
preference for somewhat lower times near ~3x 10 yrs. 

It is interesting and reassuring that the chemical evo- 
lution times from both the infrared and submillimeter 
data provide similar results. This is especially true as 
they probe such different regions of the envelope. Also, 
the nearly simultaneous agreement in time between such 
different species, with transitions arising throughout the 
envelope, and observed with a range of ground- and space- 
based instruments provides significant support to the pro- 
posed chemical, physical, and thermal structure of the en- 
velope. 



6. Conclusions 

We have constructed detailed thermal and gas-phase 
chemical models for AFGL 2591 based upon the physical 
model of van der Tak et al. ( |1999| ; |2000D . These models 
were used to probe the validity of the proposed physical 
structure, as well as study the chemical evolution of the 
source. 

In particular, we find that: 



1. Pseudo-time dependent modeling of the chemistry as 
a function of depth is a good probe of the physical and 
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Fig. 16. The quality of the model fit to the observations 
as a function of time. Here we plot the \ 2 as defined in 
the text, normalized to the maximum value for the times 
considered. The solid line shows the chi-squared for the 
submillimeter data only, while the dashed-line shows the 
results when both the submillimeter and infrared data are 
included. 



thermal structure of the source, due to the density and 
temperature dependence of the chemistry. 9. 
Care must be exercised when drawing conclusions 
from observations, as infrared absorption measure- 
ments usually better probe the warm interior of the en- 
velope while submillimeter emission data often better 
probe the cool exterior. While "average abundances" 
[x avg = iV(X)/JV (H2)] can suggest conditions differ- 10. 
ent from those actually in existence, detailed radiative 
transfer modeling are essential to give more reliable 
results (Sect. 2 & 4.1). This also underscores the im- 
portance of modeling the complete physical, thermal 
and chemical structure of the envelope when compar- 
ing to observations. 

Water and CO are stable, even at high tempera- 
tures, except for destruction by cosmic-ray-driven ion- 
molecule reactions for t > 10 5 yrs. The ice and gas 
abundances toward other sources suggest that there 
may exist a water destruction mechanism at high tem- 
peratures (> 500 K). Such a mechanism does not exist 
in our chemical model, but is worth further exploration 
(Sect. 4.1 & 4.2). 

The hydrocarbon and nitrogen chemistry is strongly 
influenced by temperature, with significant production 
possible at high (T > 400 K) temperatures. Still it 



appears that there may exist an as-yet unexplored path 
to HCN for T > 600 K (Sect. 4.2). 
The sulphur chemistry displays three significant 
regimes. In the cool exterior, over 50% of the sulphur 
exists in CS. For 100 < T(K) < 300, S0 2 contains 
a majority of the sulphur. At higher temperatures, 
atomic sulphur dominates. While consistent with ob- 
servations, these results employ a low gas-phase sul- 
phur abundance. We cannot identify the principle sul- 
phur reservoir in the cold gas (Sect. 4.3). 
Carbon dioxide may possibly be destroyed by reactions 
with molecular hydrogen in impulsive heating events 
- a reaction which requires further study. The sugges- 
tion that CO2 may be destroyed by H has problems in 
that it may destroy too much CH4 in the interior and 
produce too much atomic and molecular oxygen (Sect. 
4.4). 

The cosmic-ray ionization rate in AFGL 2591 is im- 
portant. It affects the ion abundances as well as the 
formation of HCN. These processes fix the cosmic ray 
ionization rate to better than a factor of three. Ions 
produced by cosmic rays can also effectively destroy 
H2O on timescales of > 10 5 years (Sect. 4.5). 
We produce simulated column densities of predicted 
abundant species for comparison with future observa- 
tions (Table ||). We predict, among other things, abun- 
dant atomic oxygen and nitrogen - even though oxygen 
can be shuttled to water and O2 while nitrogen is at 
most 1 % atomic with the majority in N2. The tem- 
perature variation within the source can have a large 
affect on oxygen, and hence on the rest of the chem- 
istry. (Sect. 4.6). 

While it is possible for the observed H2CO abundance 
to be reproduced by gas-phase chemistry, the same is 
not true for CH3OH. This suggests that some other 
processes are important for CH3OH, including grain- 
surface chemistry and/or the potential of a heating 
event within the envelope (Sect. 4.7). 
It is possible to use detailed chemical modeling to 
constrain the chemical age of AFGL 2591. We find 
7 x 10 3 < i(yrs) < 5 x 10 4 , with a strong preference 
for t ~ 3 x 10 4 years (Sect. 5). 

The agreement of our results with the line data lend 
significant further confirmation to the physical model 
proposed by van der Tak et al. ( 2000 ) , and the chemical 
structure and evolution proposed here. 
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